function burgers2

% code to solve, and then plot, the solution of the nonlinear BVP
%  B Euler + centered diffs
%        u_t + u*u_x = ep*u_xx     for xL < x < xR  
% where 
%          u(xL,t) = a0 ,  u(xR,t) = b0,  u(x,0) = phi(x)

%clf
% get(gcf)
%set(gcf,'Position', [1925 1095 573 199]);

% set boundary conditions
	xL=0; a0=1;
	xR=1; b0=-a0;

ep=0.01;

% parameters for calculation
nx = 400
error = 0.000001;
tmax=200
nt=200

uu=zeros(nx+2,nt);
x=linspace(xL,xR,nx+2);
t=linspace(0,tmax,nt);
for ix=1:nx+2
	uu(ix,1)=phi(x(ix));
end

plot(x,uu(:,1))
%pause

grid on
box on
%axis([-0.02 1.02 -2 2])
xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')
set(gca,'FontSize',14);


dx=x(2)-x(1);
dxx = dx*dx;
dt=t(2)-t(1);

for it=2:nt

for ix=2:nx+1
	q1=dt*uu(ix,it-1)*0.5*(uu(ix+1,it-1) - uu(ix-1,it-1))/dx;
	q2=-uu(ix,it-1);
	q(ix-1)=(q1+q2)/(dt*ep);
end;

y=uu(:,it-1);

err=1;
counter=0;
while err > error

	% calculate the coefficients of finite difference equation
	a=zeros(1,nx); c=zeros(1,nx); v=zeros(1,nx); u=zeros(1,nx);
	for j = 2:nx+1
		jj=j-1;
		z = (y(j+1) - y(j-1))/(2*dx);
		a(jj) = 2 + dxx*fy(x(j), y(j), z, q(jj), ep, dt);
		c(jj) = -1 - 0.5*dx*fz(x(j), y(j), z, q(jj), ep, dt);
		v(jj) = - 2*y(j) + y(j+1) + y(j-1) - dxx*f(x(j), y(j), z, q(jj), ep, dt);
	end;

	% Newton iteration
	v(1) = v(1)/a(1);
	u(1) = - (2 + c(1))/a(1);
	for j = 2:nx
		xl = a(j) - c(j)*u(j-1);
		v(j) = (v(j) - c(j)*v(j-1))/xl;
		u(j) = - (2 + c(j))/xl;
	end;
	vv = v(nx);
	y(nx+1) = y(nx+1) + vv;
	err = abs(vv);
	for jj = nx:-1:2
		vv = v(jj-1) - u(jj-1)*vv;
		err = max(err, abs(vv));
		y(jj) = y(jj) + vv;
	end;
	counter=counter+1;
end;

uu(:,it) = y;

newton_iterations=counter;

% plot computed solution
plot(x,y)

grid on
box on
axis([0 1 -1 1])
xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')
set(gca,'FontSize',14);

%pause
end

for ixx=1:nx+1
	if uu(ixx,it)*uu(ixx+1,it)<0
		z0 = (x(ixx)*uu(ixx+1,it) - x(ixx+1)*uu(ixx,it))/(uu(ixx+1,it) - uu(ixx,it))
		break
	end
end

%  right hand side of differential equation
function g=f(x,y,z, q, ep, dt)
g = y/(dt*ep) + q;

%  fy(x,y,z) = diff(f,y)
function g=fy(x,y,z, q, ep, dt)
g = 1/(dt*ep);

%  fz(x,y,z) = diff(f,z)
function g=fz(x,y,z, q, ep, dt)
g = 0;

%  initial condition
function g=phi(x)
x0=0.25;
a=1;
if x<x0
	g=a*(1-x/x0);
	g=a;
else
	g=-a*(x-x0)/(1-x0);
	g=-a;
end


























